
inline
gmm::gmm(const mat& in_means, const mat& in_covs, const vec& in_weights)
: means(in_means), covs(in_covs), weights(in_weights)
{

}

inline 
gmm::gmm( const char *name )
{
  load_model(name);
  
}

inline
double
gmm::likelihood(const vec& x)
{
  return std::exp( log_likelihood(x) );
}



inline
double
gmm::log_likelihood(const vec& x)
{
  // calculate the log version of probability // ( Ec. 2.20)
  double gauss_function;
  double prob;
  prob=0;
  for (uword gi = 0; gi < means.n_cols; ++gi)
  {
    gauss_function = calc_gauss_function(x,gi);
    prob = prob+weights(gi)*gauss_function;
  }
  double log_likelihood = log(prob);
  return log_likelihood;
}



inline
double
gmm::calc_gauss_function(const vec& xn, const uword gi) //Ec. 2.23
{
  double logDetCov;
  double logDen;
  double logNum;
  double logGauss;
  vec res;
  vec mu = means.col(gi);
  uword Dim = means.n_rows;
  
  double logDetCov2;
  double sign;
  log_det(logDetCov2, sign, diagmat(covs.col(gi)));
 
  //logDetCov=log(det( diagmat(covs.col(gi)) ));
  logDen = (-0.5)*(Dim*log(2*datum::pi) + logDetCov2);
  res = (xn - mu);
  logNum   = as_scalar((-0.5)*res.t()*inv( diagmat( covs.col(gi) ) )*res); //instead of cov.i(), you can do inv( diagmat( covs.col(gi) ) )
  logGauss = logNum + logDen;
  return std::exp(logGauss);
} 

inline
double
gmm::avg_log_likelihood(const mat& X)
{
  double acc = 0.0;
  
  for(uword i = 0; i < X.n_cols; ++i)
  {
    const vec tmp = X.col(i);
    
    acc += log_likelihood(tmp); // Eq. 2.19
  }
  
  if(X.n_cols > 0)
  {
    return acc / double(X.n_cols);//Ec. 2.18
  }
  else
  {
    return 0.0;
  }
}


inline 
void
gmm::save_model(const char *name)
{
  uword N_gauss = means.n_cols;
  field<mat>	gmm_model(3);
  gmm_model(0) = means;
  gmm_model(1) = covs;
  int rs = ceil(N_gauss/2);
  weights.reshape(rs,2);
  gmm_model(2) = weights;
  gmm_model.save(name);
}


inline 
void
gmm::load_model(const char *name)
{
  field<mat>	gmm_model(3);
  
  gmm_model.load(name);
  means = gmm_model(0);
  covs  = gmm_model(1);
  weights = gmm_model(2);
  weights.reshape( means.n_cols,1 );

}

